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We present a numerically feasible semiclassical (SC) method to evaluate quantum fidelity decay 
(Loschmidt echo, FD) in a classically chaotic system. It was thought that such evaluation would 
be intractable, but instead we show that a uniform SC expression not only is tractable but it gives 
remarkably accurate numerical results for the standard map in both the Fermi-golden-rule and Lya- 
punov regimes. Because it allows Monte Carlo evaluation, the uniform expression is accurate at 
times when there are 10^" semiclassical contributions. Remarkably, it also explicitly contains the 
"building blocks" of analytical theories of recent literature, and thus permits a direct test of the 
approximations made by other authors in these regimes, rather than an a posteriori comparison 
with numerical results. We explain in more detail the extended validity of the classical perturbation 
approximation (CPA) and show that within this approximation, the so-called "diagonal approxima- 
tion" is automatic and does not require ensemble averaging. 
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The question of stability of quantum motion, originally 
formulated by Peres 1], has recently attracted much in- 
terest, due to its relevance to quantum computation and 
decoherence in complex systems. Peres defined stability 
in terms of quantum fidelity Mit), the overlap at time 
t of two states, which were identical at time t — Q, but 
afterwards propagated in slightly different dynamical sys- 
tems, described by Hamiltonians and = + V , 



M{t) = |(V'|exp(iir^<) exp(-ii7°t)| 



(1) 



This quantity is also called Loschmidt echo, because it 
can be interpreted as an overlap of a state propagated 
forward for time t with and then backward for time 
t with , with the initial state. We consider to 
be strongly chaotic, although our method is not limited 
to this case. Even with this restriction, the decay of 
fidelity has a surprisingly rich behavior: Most surprising 
recently was the derivation in Ref. that for certain 
range of perturbations the decay rate is independent of 
the perturbation strength. 

The Loschmidt echo is physically realizable, for ex- 
ample in NMR spin echo experiments, where back- 
propagation under a slightly different Hamiltonian is fea- 
sible 3, 4, 5j. There arc other examples, which often go 
unnoticed. An example is neutron scattering, where the 
scattering kernel can be written as in Eq. , with a 
momentum boosted version of . Many numerical in- 
vestigations of FD have been undertaken in various sys- 
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on the strength of perturbation, there exist at least four 
qualitatively different regimes of the decay in chaotic sys- 
tems Q: As the perturbation increases, these regimes are 
perturbative (PT), Fermi-golden-rule (FGR), Lyapunov 
(L), and the strong SC regime. 

In the PT regime, in which the characteristic matrix 
element of the perturbation is smaller than the mean 



level spacing A, the decay can be described by a combi- 
nation of perturbation theory and random-matrix theory 
(RMT), and is Gaussian HQ, 



MpT {t) K, exp 



(2) 



For intermediate perturbation strengths, the decay fol- 
lows the Fermi golden rule and is exponential. 



M 



FGR 



(t) « exp(-rt/n) 



(3) 



where F = 2^]/^/ A. In Ref. [g] it is shown that this 
FGR decay is equivalent to the exponential decay derived 
semiclassically in Refs. 0,I3- other words, F — 2K/h 
where K is the classical action diffusion constant. 



K 



dt (F[r(t)]F[r(0)]), 



In the Lyapunov regime, derived in Ref. 2], FD actu- 
ally does not depend on the strength of perturbation, but 
only on the Lyapunov exponent A of the chaotic system, 



Ml (t) ~ exp {~Xt) 



(4) 



We are able to find a numerically feasible uniform 
H m SC method to evaluate FD in the FGR 
and Lyapunov regimes. As a result, we can directly test 
all approximations made in the derivation of results Q 
and ^ from Refs. 0, The method starts with a 
SC approach based on the CPA 0, Q , and ends with a 
form of initial value representation (IVR) |3^ . Is^ which 
makes the numerical calculation manageable and the SC 
approximation itself more accurate. 

Following notation of Ref. |^ , we want to find FD for 
an initial Gaussian wave packet 



i^{r;0) = (ttct^) '^^''exp 



^Po-(r-ro)-^^ 
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It is centered at Tq with dispersion a and has an average 
momentum pp. We propagate this state with a SC Van 
Vleck-Gutzwiller propagator 



K'^ir", r'; t) - ^ {27nhy''^^ C^^ exp j^S, 



n 



Here Cj = \dei Sj / dv" dr') \ is the absolute value of 
the Van Vleck determinant, Sj{r",r';t) is the action 
along the jth trajectory connecting r' with r", 

5,(r",r';i)= f dt' L[r {t') ,r (f) ,t'] 
Jo 

and Vj is the Maslov index. 

Expanding each contribution about a central trajec- 
tory 39], the overlap amplitude of the semiclassically 
propagated states becomes 



(5) 



X exp I - [(p;. - po)' + (p;-, - Po)'] (TV2ft2| 

where Sj — Sj (r, rg; t) and superscript V denotes quan- 
tities in the perturbed system. At thispoint, two crucial 
approximations are made in Refs. |3jl3' First jOnly the 
diagonal terms j = j' are considered. Ref. [3| claims 
that these are the only terms surviving the average over 
impurities in disordered systems. Below we show that 
this is not a separate approximation, but that it follows 
from the CPA and does not require any ensemble averaj 
ing. CPA, the second approximation used in Refs. 0, 
is based on an apparently hopeless assumption that the 
perturbation does not affect trajectories (i.e., Cj ~ Cj 



and w Vj) but only affects the actions, through 



Sj = - dt'V[v,{t')\. 



(6) 



Of course this assumption is wrong for individual tra- 
jectories which deviate exponentially with time. The rea- 
son why the approximation works in quantum mechanics 
is subtle: The first step to understanding why it yields 
accurate wave functions lies in the structural stability 
of the manifolds, as pointed out in Ref. 0- Assuming 
that perturbation does not cause a bifurcation and does 
not significantly change the stable manifold, the evolved 
manifolds almost exactly overlap whereas the same initial 
points deviate exponentially by sliding along the mani- 
fold 

The second step goes as follows: consider trajectories 
A{t), A^{t) under the flow iJ", iJ^, respectively. Let 
A{{)) — A^(0) be a point on the Lagrangian manifold 
supporting the wave function at i = 0. While {t) 
exponentially diverges from A[t), if the evolved manifolds 



(almost) exactly overlap, we can find a point 5(0) on the 
manifold at i = such that (t) (almost) coincides with 
A{t). Because of the exponential sensitivity to the initial 
conditions, point B{0) will be exponentially close to ^(0). 
Trajectories A{t) and B{t) remain exponentially close for 
all times, so if we use these particular trajectories to find 
i/j{t) and ip'^ (t), respectively, the CPA will be justified. 

The diagonal approximation and CPA enormously sim- 
plify expression |(SJ for the overlap amplitude: 



2\d/2 



d'^r ^ Cj exp 



(7) 



iA5,/ft- (p; -po)'aVft^ 



At this point, both Refs. 0, 13| resort to statistical ar- 
guments to obtain an analytical result. Expression Q 
for the overlap would be very difficult to implement nu- 
merically for three reasons: First, in chaotic systems 
there is an exponentially growing number of contribut- 
ing trajectories. Second, the accuracy would be com- 
promised by proliferating caustic singularities in the Van 
Vleck determinant Cj whenever dr/dp'^ — 0. Finally, for 
each trajectory we would have to perform a computation- 
ally expensive root-search to find initial p'j that satisfies 
r(ro,Pj-,i) = r. However, there exists a beautiful and 
simple way to eliminate the exponential number of con- 
tributions, caustic singularities, and the root search, all 
at the same time. All three problems can be solved if we 
evaluate the overlap Q in the initial momentum instead 
of final position representation. Exactly one point on the 
evolved manifold corresponds to each initial momentum, 
so no summation is necessary. The new "Van Vleck de- 
terminant" is exactly 1, so there will be no Maslov indices 
either. With all these simplifications, the SC evaluation 
becomes tractable; in principle, it yields the same result 
that an arduous evaluation of iQ would: 



Oit) = (aVvrft^)'^/^ / 



(8) 



exp iAS[r{rn,p',t),rn,t]/h-{p' -pof<T^/h^ . 



The only assumption required to derive © is the validity 
of CPA, in the extended sense described above. Ensem- 
ble averaging used in Ref. is unnecessary: result © 
works for pure states. Expression ||SJ| is a special form 
of IVR jSy, [33 • In general, IVR avoids the singularities 
and the root search, but at a cost of replacing a sum 
over classically allowed paths by an integral over all ini- 
tial momenta. In our case, it is even better, since we also 
eliminated the integral over final position r. We remark 
that ^ can also be obtained by changing the integration 
variable in |(7J) from final r to initial p', but our deriva- 
tion avoids the intermediate step {Tj) that requires mak- 
ing diagonal approximation in (0. We note the unique 
property of IVR: in this representation, FD is only due to 
dephasing. In other representations, the decay can also 
have a component due to the decay of classical overlaps. 
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We chose to test our method on the standard map used 
in Ref. i, 

1j+i=1j+P3 (modi), 
k 

Pj+i =P3 - TT sin (27r(jj+i) (mod 1). 
zvr 

Perturbation is effected by replacing the parameter k by 
k + e. Choice of an n-dimensional Hilbert space for the 
quantized map fixes the effective Planck constant to be 
h = (27rn)~^. We note that results of exact quantum 
and SC computations which we present below are for 
initial position eigenstate with qo = 0.5 rather than a 
wavepacket. 

In previous numerical experiments analytical predic- 
tions of Gaussian or exponential decay have been com- 

Pired to an exact quantum calculation: see, e.g., Refs. 
While we also have an exact quantum 
benchmark (FFT) with which to compare the expres- 
sions for various regimes, we reiterate that it would be 
hard from a mere comparison of final results for AI (t) to 
determine the source of errors. We proceed by discussing 
how the uniform method helps to analyze various regimes 
of decay. In the PT regime (see Fig. ^ , we do not expect 
any SC approach to work very well except for short times 
(much shorter than the Heisenberg time iff = h/ A). The 
RMT analytical result MpT from Ref. [3 gives an ex- 
cellent agreement in this case. The inset shows, how- 
ever, that before the Gaussian decay MpT sets in at the 
Heisenberg time, the uniform expression follows M^xact 
much better. 
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FIG. 1: Fidelity in the perturbative regime (fc = 18, A « 2.21, 
e — 10"'', tH ~ ri = 350). Inset: detail for short times. 

As the perturbation e increases, we enter the regimes 
with exponential decay of fidelity. If the perturbation is 
strong quantum mechanically, but does not significantly 
change the stable manifold, CPA may be used. Even 
within the CPA, there are two types of decay, discussed 
already in Ref. 0]. First, there is decay related to de- 
phasing of trajectories with uncorrelated actions. Sec- 
ond, there is decay related to dephasing of very near 
trajectories with correlated actions. For smaller per- 
turbations, the first type of decay is slower and domi- 
nates the behavior of fidelity: this happens in the FGR 



regime. For larger perturbations, dephasing of uncor- 
related trajectories is so fast that the quantum overlap 
is determined by the fraction of near trajectories that 
have remained in phase. This is the case in the Lya- 
punov regime. Transition from the PT to FGR regime 
occurs for w 32Tr^n^^[l + 2J2{k)]~^ 7] when most of 
the overlap has decayed before Heisenberg time. Tran- 
sition from the FGR to Lyapunov regime occurs for 

« Sn'^Xn^^ll + 2J2(fc)]"i when the FGR decay rate is 
larger than A. 

Using Eq. 0, fidelity can be written as a weighted 
average of terms exp [i{AS' —AS")/h], 



Mumf (^)-(^ 



d'^p'J d'^p" exp 
exp {- [(p' - Po)' + (p" - Po)'] a^h^} 



-{AS'^AS") 
h 



(9) 



where AS"' corresponds to a trajectory with initial mo- 
mentum p". Assuming the averaging window (i.e. the 
momentum width of the wavepacket) is large enough, we 
can make the replacement 

exp [i{AS'-AS")/h] « (exp [i{AS' - AS") / h]) (10) 

in Eq. where averaging is over all initial momenta p' , 
p" . In the FGR regime where dephasing is determined 
by uncorrelated trajectories, a further simplification 

(exp [i{AS'-AS")/h]) « (e'^s'/'') (e-'^^"''') (11) 

is possible. Due to the central limit theorem, in chaotic 
systems distribution of AS approaches a Gaussian and 

(exp (iAS/h)) = exp [i {AS) /h - o-is/2?i^] (12) 

where a'^g = 2Kt is the action variance at time t. Apply- 
ing approximations (|l()|l . Ijllj) . p2|l in Eq. 0, confirms 
Eq. © for the FGR decay (a|7[. Fig. H shows FD in the 
FGR regime. In the inset, histogram of action differences 
is compared with a Gaussian fit, confirming assumption 
(|12ll . It is apparent that M^nif matches M^xact better 
than the MpcR since M^nif takes into account the pre- 
cise initial conditions without the averaging assumption 
(|10|l and since MpGR uses an analytic result for K , which 
is only approximate . Careful inspection of the short 
time regime (not shown) reveals that M^nif agrees with 
Mexact, since unlike Mpqr, Munif does not depend on 
the central limit theorem which guarantees the Gaussian 
assumption (|12|l at later times. Finally, we would like 
to point out that the uniform expression is very accu- 
rate at time t « 120 when there are approximately lO'''^ 
semiclassical contributions in the sum Q! 

In the Lyapunov regime, FD is determined by dephas- 
ing of near trajectories with correlated actions |3, in- 
validating simplification (|11() . Now the action difference 
AS"— AS"' depends on the initial momenta p', p" . Using 
reasoning similar to Ref. 2] or statistical arguments for a 
random walk with an exponentially increasing time step 
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FIG. 2: Fidelity in the FGR regime {k = 18, A « 2.21, 
e = 5 X 10"'', n = 3500). Horizontal dashed line ("ergodic") 
is the limit of FD due to the finite size of Hilbert space. Inset: 
Histogram of action differences compared to a Gaussian fit. 



|4(3|. it can be shown that the action difference is also 
Gaussian distributed, with zero average and variance 

([AS ip') - AS {p")f) « {D/2X)e'^' (p' - p"f , (13) 

poo 

D^2 dt {V {Q)V' {t)) . 
Jo 

We can therefore make the replacement 

D 



^exp 



:(AS"-A5" 



cxp 



e^'^'ip'-p"? 



in Eq. © and ^ to find 
Ml it) « (l + e2^*2?/2Aa2)"'/' 



-At 



confirming Eq. For the precise definition of A, see 
Ref. ^3 (one has to be careful about the averaging pro- 
cess). Fig. 121 displays M {t) in the Lyapunov regime. It 
shows that while Ml gives an accurate average decay 
only for At 3> 1, M^nif correctly follows the behavior of 
Mexact even for short times t A^^. The inset shows 
the variance of AS* {p') — AS [p") as a function of p' —p" 
at a fixed time and justifies the assumption made in Ref. 
in derivation of perturbation independent decay: for 
near trajectories, the variance grows quadratically with 
p' ~ p" (fitted line gives an exponent 2.003), in accor- 
dance with Eq. (|13(l . while for distant trajectories, in 
accordance with the derivation of the FGR regime, the 
variance is independent of p' — p", 



[AS [p') - AS {p")f) = 2als = ^Kt. (14) 



The time dependence of [AS (p') - AS ip")fj for fixed 

p' — p" is shown in Fig. 01 Part a) shows that for short 
times when trajectories are still correlated, this depen- 
dence is exponential, in agreement with Eq. (|13|l . Part 



b) shows that for longer times, when correlation is lost, 
the dependence is linear, as expected from Eq. (|14|) . 
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FIG. 3: Fidelity in the Lyapunov regime (fc = 7, A ~ 1.28, 
e = 5 X 10"*, n = 10^). Meaning of lines same as in Fig. |21 
Inset: Variance of AS (p") — AS (p') as a function of p" — p' 
at time t = 7. Dots are numerically calculated, dashed line 
is the horizontal asymptote 2a\s^ solid line is a linear fit for 
small p" — p", in agreement with Eq. 11311 . 




FIG. 4: Variance of AS (p") - AS (p) as a function of t for 
p" — p' = lO"^'^: a) exponential dependence for short times, 
b) linear dependence for long times. 



To conclude, we have explicitly evaluated SC expres- 
sions which were thought to be intractable numerically, 
yielding remarkably accurate results for FD in the FGR 
and Lyapunov regimes. We provided a more detailed ex- 
planation why CPA works and employed our method to 
test other approximations used in Refs. 0, ■ 
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